% Calculate and plot similarity matrices
% Plot basis functions and energy in modes
clear all; close all; clc;

species = {'Abalone';'Lpictus';'Lvariegatus';'Spurpuratus'; ...
    'Arbacia';'Ciona';'Zebrafish';'Bull';'Human-low';'Human-high';};

use_optimization = 1; % 1 - use optimization, 0 - don't use optimization
dat = struct([]);
recalcAll = false;

date = {'200617B'; '200617B'; '200617B' ; '200617B'; ...
    '200617B'; '200617B'; '200617B' ; '200617B';...
    '200617B'; '200617B'};
savedate = '200629-KD';

colors{1} = [0 1 1];
colors{2} = [1 0 0];
colors{3} = [1 0 0];
colors{4} = [1 0 0];
colors{5} = [1 0 0];

colors{6} = [1 0 1];
colors{7} = [0 1 0];
colors{8} = [0 0 1];
colors{9} = [0 0 1];
colors{10} = [0 0 1];

actualSpecies{1} = 'H. rufescens';
actualSpecies{2} = 'L. pictus';
actualSpecies{3} = 'L. variegatus';
actualSpecies{4} = 'S. purpuratus';
actualSpecies{5} = 'A. punctulata';
actualSpecies{6} = 'C. intestinalis';
actualSpecies{7} = 'D. rerio';
actualSpecies{8} = 'B. taurus';
actualSpecies{9} = 'H. sapiens';
actualSpecies{10} = 'H. sapiens (viscous)';

symb = 'ovs^oooovp';

for x = 1 : length(species)
    
    filename = strcat(species{x},['_' date{x} '.mat']);
    load(filename)
    dat(x).Uall = U_temp;
    dat(x).Sall = S_temp;
    %if exist('procdata')
    %    dat(x).U = double(procdata(1).U);
    %    dat(x).S = double(procdata(1).S);
    %    clear procdata allruns U_temp thresh kfolds
    %else
    U_temp = cat(3,U_temp{:});
    
    stdU{x} = std(U_temp,[],3);
    U_temp = mean(U_temp,3);
    sEnd = size(S_temp{1},1);
    for i =1:length(S_temp)
        try
            S_temp{i} = S_temp{i}(1:sEnd,1:sEnd);
        catch
            sz = size(S_temp{i});
            S_temp{i} = [S_temp{i}(1:sEnd,1:sz(2)), zeros(sEnd,sEnd-sz(2))];
        end
        svi = diag(S_temp{i});
        dat(x).Sall{i} = double(cumsum(svi.^2) ./ sum(svi.^2));
    end
    S_temp = cat(3,S_temp{:});
    S_temp = mean(S_temp,3);
    sv = diag(S_temp);
    dat(x).U = double(U_temp);
    
    dat(x).S = double(cumsum(sv.^2) ./ sum(sv.^2));
    
    
    clear A_overthresh allruns U_temp thresh kfolds S_temp
    %end
end

%% Kinematic Error on all individual curves
%wMval = zeros(10,10,sum(1:49));
wMval = nan(10,10,100);

s = RandStream('mlfg6331_64');
for p = 1:100
    rands100(p,:) = datasample(s,1:50,2,'Replace',false);
end

for i = length(dat):-1:1
    for j = length(dat):-1:1
        counter = 0;
%         try
%             loadfile = load('200617-EL2-KD-Final_ALLwMerror.mat');
%             wMval(~isnan(loadfile.wMval)) = loadfile.wMval(~isnan(loadfile.wMval));
%         end
        try
            loadfile = load([savedate '_ALLwMerror.mat']);
            wMval(~isnan(loadfile.wMval)) = loadfile.wMval(~isnan(loadfile.wMval));
        end
        %for ind1 = 1:50
        tic
        for p=1:50
            try
                ind2 = rands100(p,1);
                ind1 = rands100(p,2);
                %for ind2 = (ind1+1):50
                counter = counter+1;
                U_{1} = dat(i).Uall{ind1}; U_{2} = dat(j).Uall{ind2};
                %S_{1} = dat(i).Sall{ind1}; S_{2} = dat(j).Sall{ind2};
                S_{1} = dat(i).Sall{ind1}; S_{2} = dat(j).Sall{ind2};
                if isnan(wMval(i,j,counter)) || recalcAll
                    wMval(i,j,counter) = calculateWMval(U_, S_);
                end
                
                while wMval(i,j,counter) > 1 || wMval(i,j,counter) < 0
                    [ind1, ind2] = datasample(s,1:50,2,'Replace',false);
                    U_{1} = dat(i).Uall{ind1}; U_{2} = dat(j).Uall{ind2};
                    %S_{1} = dat(i).Sall{ind1}; S_{2} = dat(j).Sall{ind2};
                    S_{1} = dat(i).Sall{ind1}; S_{2} = dat(j).Sall{ind2};
                    wMval(i,j,counter) = calculateWMval(U_, S_);
                end
                %end
            end
            tempvals = wMval(i,j,:);
            disp(['Mean through ' num2str(counter) ' of ' num2str(100) ': ' num2str(nanmedian(tempvals(tempvals~=0)))]);
            disp(['Std through ' num2str(counter) ' of ' num2str(100) ': ' num2str(nanstd(tempvals(tempvals~=0)))]);
            
            if mod(p,10)==0
                save([savedate '_ALLwMerror.mat'],'wMval');
                figure(1); imagesc(nanmedian(wMval,3)); axis image; colormap hot; colorbar;
                figure(2); imagesc(100-sum(isnan(wMval),3)); axis image; colormap gray; colorbar;
            end
            
        end
        toc
        
    end
end

%medWM = nanmedian(wMval,3);
% KD_ = medWM; 
% KD_(2,:) = medWM(4,:); KD_(:,2) = medWM(:,4);
% KD_(4,:) = medWM(2,:); KD_(:,4) = medWM(:,2);
% KD_(2,4) = medWM(4,2); KD_(4,2) = medWM(2,4);
% KD_(2,2) = medWM(4,4); KD_(4,4) = medWM(2,2);
% medWM = KD_;

medWM = nanmedian(wMval,3);
wMij = mean(cat(3,medWM,medWM'),3);
figure; imagesc(wMij); axis image; colormap hot; colorbar;
set(gca,'xtick',[1:10],'xticklabel',actualSpecies,'ytick',[1:10],'yticklabel',actualSpecies);
xtickangle(90);

figure; imagesc(wMij(1:9,1:9)); axis image; colormap hot; colorbar;

for i=1:10
    for k=1:10
    bpWMt{i}{k} = squeeze(wMval(i,k,:)); 
    end
end
figure; violin(bpWMt{1},[],'labels',actualSpecies,'Color',[0.0 0.5 0])
